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In this paper, we study transport features of a one-dimensional beam-plasma system in 
the presence of multiple resonances. As a model description of the general problem of a 
warm energetic particle beam, we assume n cold supra-thermal beams and investigate the 
self-consistent evolution in the presence of the complete spectrum of nearly degenerate 
Langmuir modes. A qualitative transport estimation is obtained by computing the 
Lagrangian Coherent Structures of the system on given temporal scales. This leads to the 
splitting of the phase space into regions where the local transport processes are relatively 
faster. The general theoretical framework is applied to the case of the nonlinear dynamics 
of two cold beams, for which numerical simulation results are illustrated and analyzed. 


1. Introduction 

Nonlinear interplay of energetic particles (EPs) with Alfvenic fluctuations, such as 


Alfven eigenmodes (AEs), EP modes and drift Alfven waves (Chen & Zonca 2015 


2007, 2013 Fasoli et al. 20071, and their consequences for fluctuation induced cress¬ 


Heidbrink||2008 Breizman & Sharapov 

2011| Lauber 2013 

Vlad et a/.|2013 Wang et al. 

2mij 

Schneller et al. 2013 2012 Gore 

enkov et al. 2014 

Zonca et a7|2015a|6| 

Pinches 


et al. 


20151. The relevance of the beam-plasma instability problem as paradigm for EP 
transport due to AEs was proposed in ( Berk fc Breizman||1990a|()|c ). In fact, the bump- 
on-tail paradigm can mimic the nonlinear interaction of supra-thermal particles with 


Alfvenic fluctuations, provided the system is sufficiently close to marginal stability (Chen 
fc Zonca|2007 2013[ Berk &: Breizman|1990c Breizman &: Sharapov|20fT ). In particular, 
it is possible to determine a correspondence between velocity space transport in the one¬ 
dimensional (ID) beam-plasma instability problem and radial transport of fast ions in 
the presence of resonant AEs, at least as far as their spatial structure does not play 
a significant role. More precisely, such correspondence has rigorous foundation as long 
as the nonlinear resonant particle displacement is small compared to the characteristic 
perpendicular fluctuation scale (Chen & Zonca 2015 Zonca et a7||2015a|6 Wang et al. 
2012 Briguglio et aL|[2M4 |. 

The beam-plasma instability was analyzed in the pioneering works of O’Neil and 
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collaborators (O’Neil et al. 1971 O’Neil & Malmberg 1968) [see also (Shapiro 1963 
Onishchenko et al. |1970[ [Shapiro fc Shevchenko|1971[|Thoinpson|197H [Levin et a/.jl972| 


Matsiborko et al. |1972[ )]. It was analyzed for a nearly monochromatic fast electron beam 
injected into a ID plasma, which, in turn, is treated as a linear dielectric medium sup¬ 


porting longitudinal Langmuir waves (Lifshitz & Pitaevskii 1976). It was demonstrated 


that Langmuir fluctuations can be resonantly excited by EPs and become unstable. 


The original approach of ( 

O’Neil et al. 

1971 

) can be cast in a rigorous Hamiltonian 

formulation ( 

Mynick & Kaufman 

1978 

Tennyson et al. 

1994 Antoni et al. 

1998 

); and 


the resulting scheme is formally equivalent to that of the Free Electron Laser dynamics 


Antoniazzi et al. 

2008 

). An alternative approach was given in ( 

Al’Tshul’ & Karpman 


1966), where the nonlinear Vlasov-Poisson initial value problem was solved for the case 


of nearly degenerate Langmuir and beam modes and an analogue of the Dyson equation 
was derived [see ( Chen fc Zonca|2015 Zonca et aZ.|2015a|l) ) for an in-depth discussion of 
this approach]. An example of the beam-plasma instability in the presence of more than 
a single resonant mode was analyzed by (Zaslavsky et al. 2008), emphasizing the link 


with the so-called sideband instability, and addressing a simple Hamiltonian formulation 
of the problem. Meanwhile, ( Evstatiev et a7||2003 ) also include the effect of a magnetic 
field in the multi-wave model, but the analysis refers to a single longitudinal resonant 
electrostatic mode. Pioneering works on these thematics can be found in ( Deneef|[l975 ) 
and in (Lacina et a7|1976), where the ID beam-plasma interaction with higher harmonics 


has been studied. 

The beam-plasma problem is characterized by typical velocity and temporal scales. 
Nonlinear saturation of the beam-plasma instability occurs over a ^ characteristic 
time, with the linear instability growth rate ( O’Neil &: Malmberg||1968 ); while beam 
heating typically takes place at a rate lob oc i.e., the wave-particle 

trapping frequency in the fluctuating scalar potential Lp with wave number k (here, rrif, and 
e denote the mass and charge of electrons, respectively) and involves resonant particles 
in an interval Av ~ aJs/k ~ {u!b/uJp)vb about the velocity resonance condition v = vb 
(where vb the EP drifting beam velocity). The saturation level of the beam-plasma 
instability can be estimated as wb ~ 37b ( [O’Neil et a7|197T Levin et a7|1972 ), while the 
flattening of the beam distribution function can take much longer times than ^ Wg 


and 


involves the formation of structures in the wave-particle trapping region (Tennyson et al. 


1994). It has been recently shown (Carlevaro et al. 2014) that the nonlinear evolution 


of the beam-plasma system is characterized by quasi-stationary states which can be 


accurately predicted by the maximum entropy principle proposed by Lynden-Bell (Antoni 

erZ][T9^ . 


In the present work, we formulate the general problem of n cold beams self-consistently 
evolving in the presence of m ^ n unstable modes degenerate with the Langmuir 
potentials at the plasma frequency uip, i.e., ui = ujp, which ensures that the dielectric 
function of the cold background plasma is nearly vanishing. Among these waves, we 
can select n modes which are resonantly driven by the beams [obeying to the O’Neil 
selection rule k = uJpjvB (O’Neil & Malmberg 1968 O’Neil et a7)1971 )] and m — n quasi¬ 
degenerate ones associated with nonlinear wave-particle and/or mode-mode coupling in 
the case of (almost) overlapping resonances. It is worth noting that the dominant quasi¬ 
degenerate modes can also be dealt with as slowly evolving Langmuir waves when the 
spectral density is sufficiently large. In particular, we will show in Sec. [^ that nonlinear 
excitation of quasi-degenerate sidebands can be important for the case of two interacting 
beams. 

The physical motivation for the proposed model relies on the fact that, due to the 
intrinsic beam heating involved in the nonlinear beam-plasma dynamics (Shapiro|1963), 
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transport is not qualitatively different for cold and warm beams as long as the initial beam 
width in velocity space is smaller than Av. Thus, transport in the presence of overlapping 
resonances may be studied as the evolution of initially cold beams as a consequence of the 
corresponding nearly degenerate beam (Langmuir) instabilities. This statement justifies 
the representation of the bump-on-tail problem in terms of n overlapping resonances and 
has the valuable feature of maintaining direct control on spectral density and intensity, 
which are the properties that ultimately determine EP transport. 

The proposed theoretical model can properly describe the case of isolated resonances 
and m = n uncoupled nonlinear oscillators (dominated by wave-particle trapping), as 
well as the “quasi-linear” diffusion regime typical of a “broad fluctuation spectrum” 
( Al’Tshul’ fc Karpman|[l966 ). In fact, the nonlinear equation system we propose, in the 
limit of very large n, is expected to be physically isomorphic to quasi-linear models of EP 
interaction with AEs. In this respect, it is important to mention that a line-broadened 
quasi-linear approach has been proposed by (Berk et al. 1995a 1996a) for computing 


EP transport by means of a diffusion equation, which could address not only overlapping 
modes but also the broadening of the resonant spectrum for isolated instabilities in the 
case of multiple AEs. This model has been extended and compared with experimental 


observations (Ghantous et al. 2012) and with numerical solutions of the bump-on-tail 
paradigm (Ghantous et al. 20I4[). Contrary to th is approach, which cannot account for 
frequency chirping phenomena (Ghantous et al. 2014), the present n-beam model aims 


at characterizing particle transport on time scale comparable to the instability growth- 


rate and, thus, it can address a broad range of underlying complex phenomena (Chen & 
Zonca|[^015 Zonca et aL||2015a|6 ). 


As an illustration of the onset of particle transport phenomena in the phase space, the 
simulation discussed in this work are specialized to two cold beams and corresponding 
instabilities. Clearly, when the two beams have sufficiently different initial velocities, we 
observe the growth and saturation of two independent modes and transport processes 
are practically suppressed. Meanwhile, when the two velocities become closer, a regime 
exists when the slower beam signihcantly interacts with the resonant mode excited by 
the faster beam, and an 0(1) transfer of energy from the latter mode to the former 
takes place. Einally, when the velocity ratio exceeds a critical threshold consistent with 
the Chirikov criterion ( Levin et aL||1972 Chirikov]|1979 ), significant phase-space mixing 
takes place between the two beam populations. 

As we are interested in the details of the the mixing process between supra-thermal 
particles belonging to different beams in the presence of unstable fluctuations, we will 
present an analysis using a Lagrangian approach in order to quantitatively define the 
transport barriers of the system by means of the so-called Lagrangian Coherent Structures 
(LCSs), introduced for the first time by ( Peacock fc Haller|2013| ) and defined in terms of 
the features of the Finite Time Lyapunov Exponent (FTLE) fields ( Falessi et a/.||20f5 ). 
For practical reasons, particle transport is described as the evolution of single passive 
tracers (test particles) in a time-dependent potential generated at each instant by the self- 
consistent charge distribution. This approach, which is equivalent to the original problem, 
leads to important conceptual simplifications if we are interested in studying the onset 
of an instability that naturally selects only some Fourier components. The advection of 


passive tracers in a non-periodic flow has been studied by (Shadden et al. 2005 Haller 


2011), leading to the definition of LCSs. These structures, which are ID curves in our 


system, are hnite-time transport barriers for the tracers and, therefore, play the role of 
the last periodic surface of the Poincare map introduced in ( Shadden et aL|[2005 Haller 
2011). Thus, they divide tracers into two classes with a qualitatively different evolution. 


This way of characterizing transport gives a criterion for the transition to local chaotic 
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behavior, i.e., for assessing when nonlinear dynamics causes neighboring cold beams to 
merge. 

The paper is organized as follows. In Sec[^ the nonlinear treatment of the beam- 
plasma system of O’Neil is generalized in the presence of n distinct beams and m nearly 
degenerate Langmuir modes. In Sec|^ the theoretical model is then specialized for the 
two beam case and, in the first subsection, for the corresponding two most unstable 
modes. Numerical simulations are presented first by analyzing temporal evolutions of 
scalar potentials and then by discussing the phase-space dynamics. Mixing phenomena 
between different populations are qualitatively analyzed as a function of the ratio be¬ 
tween initial beam velocities. An analysis of quasi-degenerate modes including nonlinear 
sideband excitation is discussed in the second subsection. In Sec|^ transport processes 
are qualitatively studied by means of the LCS technique. Phase-space mixing is found to 
be enhanced at the intersections between LCSs belonging to different beam populations. 
This analysis allows us to define a critical initial velocity ratio of beams as the smallest 
value leading to beam mixing. Concluding remarks are given in Secj^ 


2. Nonlinear evolution of n resonant beams 

In the original stability analysis of the beam-plasma system ( O’Neil fc Malmberg|1968 ), 
it was shown how fluctuations can be resonantly excited by an injected supra-thermal 
beam (of velocity vb)- Introducing the plasma frequency Wp = ^Fnnpe^ jnie {rip being 
the plasma density), two types of modes were shown to exist as solutions of the dispersion 
relation: (i) the stable Langmuir wave with oj = Wp, modified by the presence of EPs; 
and (ii) the beam mode, with both stable and unstable branches at frequency ui ~ kvB, 
and maximum growth rate when uj ujp (here, k represents the wave number of the 
unstable longitudinal mode). Thus, most unstable beam modes are excited when they 
are “nearly degenerate” with Langmuir waves, i.e., k c s ujp/vB- Mean while, the nonlinear 
evolution takes place in two stages ( O’Neil et al.|197T Shapiro|1963 ): first, beam-plasma 
interactions heat the beam as EPs begin sloshing back and forth in the potential well of 
the wave; second, the beam distribution function is eventually flattened in velocity space 
by nonlinear interactions. Thus, the process leading to saturation of the beam-plasma 
instability is due to wave trapping of resonant particles, whose transport in velocity space 
is connected with energy exchange between fast electrons and the electric field. 

In this Section, we derive the governing equations for the general system of n cold 


beams and n self-consistently coupled modes (O’Neil et al. 19711. It is important 


to stress that all considered m modes are assumed to be nearly degenerate Langmuir 
waves with close to vanishing dielectric function, while resonance conditions are not 
imposed a priori. The background plasma is a cold homogeneous linear dielectric medium 
with constant particle density rip. The dynamics of fast electrons is determined by the 
Newton’s law. Equations of motion are closed by the Poisson equation which provides 
the time evolution for each of the nearly degenerate Langmuir waves considered in the 
model. 

In a ID model, the motion along the x direction of period L is labeled by xi, X 2 , •.., 
Xn for each beam, respectively. Decomposing the single beams in A^i, N 2 , ..., iV„ charge 
sheets located at xu, X 2 i, ■.., Xni, the system can be discretized using the following total 
charge density 

PB{x,t) = - Xu) + ... + -^^S{x - ; 


( 2 . 1 ) 
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where ub is the total supra-thermal particle number density while aius, crnUs are 
the number densities of each beam. The Langmuir wave scalar potential can be expressed 
using the Fourier representation as 

m 

^{x,t) = + C.C.) , (2.2) 

i=i 

where (pj {kj , t) are assumed to be slowly varying fields and we have explicitly indicated 
the time dependence relative to the rapid oscillations associated to the plasma frequency 
(the c.c. denotes the complex conjugate). Thus, the (/?j’s represent the effective nonlinear 
evolution of the potentials. 

The trajectories of charge sheets are derived from the motion equation x = eVtp/rrie 
as 

ra 

xi, = - V + c.c. , (2.3a) 

m 

i=i 


'^ni — ^ 

m ^ 

i=i 




(2.36) 


Meanwhile, the evolution of the fields (pj is determined by the Poisson’s equation written 
in the Fourier space for fixed kj as 

= AnpBjikjN) , (2.4) 

where Cp is the dielectric function (we remind that for a cold plasma Cp = 1 — Wp/w^) 
and pBj(kj,t) is the Fourier transformed density 


1 


PBj(kj,t) = 7 / dx pBix,t)e 


—iki 


= -euB 


Ni 


, — ikjXii 


+ ... + — 

N 

dy n. 


Nn 


(2.5) 

For a cold plasma and j = 1, 2, ..., m nearly degenerate Langmuir waves, the dielectric 
function is nearly vanishing and can thus be expanded near (wp, kj) as 

€p{uj{t),kj) = €p{ujp,kj){dcp/duj) {ujj - ujp) = e'{ujj - ujp) , ( 2 . 6 ) 

ujp,kj 

where the time-dependent ujj account for the nonlinear frequency shift with respect to 
ujp. Without loss of generality, the fields (pj can be now expressed as 

(2.7) 


p>j{t) oc exp(—i / dt{ijjj{t) — ojp)) , 

Jo 


thus obtaining ipj = —iipJj — ujp)ipj; and Eq.(2.4) finally writes 


ipj = 




A 

. 47renB / ui 


—ikjXii-\-iLOpt 


dy n. 


Nr, 


, — ikjX\i-\-iuj. 




( 2 . 8 ) 


The equations of motion governing the beam-plasma system are Eqs.(2.3) and (2.8). 

Let us now analyze the effective nonlinear shift ( of charge sheets due to the beam-wave 
interaction by decomposing the motion as 


X 


ni — 


Xu = VitCu{t) , X2^ = V2tC2i{t) , 


Vnt + Cni{t)- ( 2 . 9 ) 
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Performing this coordinate transformation and introducing the scaled variables (we recall 
that, for a cold plasma, e' = 2/wp) 

kj = 2n£jlL , = 2TTCni/L , r = tUpf] , (j)j = [ipj efc|]/[TO^^a;p] , (2.10) 

with fj = (jiBI2.np)^l'^ and the frequency mismatches 

Plj = “ WpJ/WpT? , ... , Pnj = [kjVn “ Wp]/Wp77 , (2-11) 

the set of equations describing the effective nonlinear dynamics of n cold beams inter¬ 
acting with m Langmuir modes reads now (the prime denotes the r derivative) 

m 

Cu = * E ^7' + c.c. , (2.12a) 

i=i 


i=i 

Ni 


+ c.c. 




-iijCli-ihjT , 


-‘■^71. 


Nrr 




)■ 


(2.126) 

(2.12c) 


The frequency mismatch terms are crucial (as in the single wave model) for defining which 
mode is growing fastest. The leading modes correspond to the ones which minimize the /?’s 
implying that the most unstable modes associated to a given beam must obey the O’Neil 
selection rule kj = Vjjojp. We stress how this dynamical system is isomorphic with respect 
to the model commonly adopted in a class of studies of wave-particle interactions and 
rigorously derived in ( Elskens fc Escande||2003|) (see a lso ( Antoni et al.||1998[ Antoniazzi 
et al. 2006 Firpo & Elskens 2000 Farina & Pozzoli 2004 Farina et al. 19941). For an 
example of the mapping between the beam-plasma instability and the Free-Electron Laser 
model see ( Antoniazzi et aZ)]|2005 ). 

The conservation laws of the system above can be derived by extension of the method 
outlined by (O’Neil et a/.||197T). The total momentum conservation reads now 


E 

i=i 




E *^7 + - + 7^ E - I^M , 


AT, 






(2.13) 


with constant /Cm- As already stated in Eq.(2.7), the potentials obey to the equation 
ipj = — ujp)Lpj, which can be simply rewritten as (p'j = —if2j(f>j, where fij = 

{uij — ujp)/(copfj) = Sujj/ujpfj. Using these expressions, the total energy conservation can 
be cast in the following form 


l<^il 


y-IU 

fi 






1 

2 LlVi 


i=i 

Ni 

, r,_.. . n2 




Ntz 


E(‘^7 + 27rwi/L7?a;p)2 -f 


... -f 






= /c 


E , 


(2.14) 


with constant /Ce- These conservation laws are directly connected with the Hamiltonian 
representation of the system. Thus the present approach based on the motion of IV 
particles in the self-consistent Langmuir potentials is more general and reduces to the 
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Vlasov-Poisson description in the mean field limit of the Hamiltonian formulation treated 
a la Klimontovich (Klimontovich 19671. The Vlasov-Maxwell equations are exactly 
recovered taking the continuum or mush limit, i.e., artificially increasing the number of 
particles and keeping the charge density constant. A discussion of this paradigm in view 
of the Lynden-Bell treatment of the distribution function for the beam-plasma model is 
presented in ( Carlevaro et a/.|[2014 Antoniazzi et aL|2008 |. The significant advantage in 
adopting the present discrete model relies on the simplicity of the governing equations for 
particle and field self-consistent evolution. Furthermore, the present scheme allows direct 
control of the spectral density {i.e., the separation in velocity space of the beams) and 
on the spectral intensity {i.e., the values of the cr’s) of the initial beam-plasma system. 

As to the comparison of the present analysis with respect to the Hamiltonian approach 


presented in (Elskens &: Escande 2003 Escande &: Elskens 20081 [see also (Zaslavsky et al. 
2008 Krafft & Volokitin 2014 Volokitin & Krafft 20121 and references therein], it is worth 
noting that they are closely related except for their analysis of an assigned warm beam 


in contrast to our approach, based on n cold beams. Due to the analysis in (Levin et al. 


19721, the intrinsic beam heating by nonlinear wave-particle interactions extends up to a 
velocity width of the order Av ~ ws/A: ~ ~ 3V5fju}p/2k. When multiple resonant 

beams overlap, the formation of a continuous particle distribution function is simply due 
to the velocity spread associated to the beam heating. Thus, addressing the hot beam 
dynamics by n cold overlapping resonant beams is a reasonable and reliable simplification 
of the nonlinear beam plasma dynamics. 

Let us now derive the linear dispersion relation [see ( O’Neil fc Malmberg|1968' | for the 
analysis of the linear beam-plasma instability]. As first step, we assume J7j independent 


of r, leading to the expression cfj = cjNe 


0 


The equation of motion can then be 


solved using linear perturbation theory by expanding (fai = Cai + Ci where 

the subscript denotes the “a” beam. The first-order integration yields 


Cai — Cai + Ca° T - » 

i=i 


(■J Vi exp[i£j-Cai + + iPajT] 




-b c.c. 


(2.15) 


The Poisson equation (2.12c) can be integrated by substituting the expression above and 


by considering uniform initial conditions for the nonlinear shift such that = 0. 

Expanding the exponential in the small (f>j limit and neglecting 0{(j)^), we finally obtain 


= —iilj(f>j = 




— i 


Na 


Na 

E 

2 = 1 




(2.16) 


According to Eq.(2.9), we consider = 0; thus, the equation above reduces to the 


following dispersion relation 


[2j — 




0-2 


_i_ -I__ (2 lyi 

(12, - /3i,)2 ' 1 1 

where the presence of the j-index represents the corresponding mode for which the 
relation is evaluated; and we recall that /3,-, = 0. 

To conclude this Section, since the proposed model (as discussed above) incorporates 
the bump-on-tail paradigm, let us now briefly discuss some features and issues of this 
problem. The bump-on-tail paradigm and its generalization in the presence of dissipation 


and collisions (Breizman et al. 

1993 Berk et al. 

1995b 

19966 1 

997 I 

Breizman et al. 

1997 

Wong & Berk 

1998 

) was extensively studiec 

in ( 

Lesur et al 

2009 

2010 

Lesur & 


Idomura 2012), along with an in-depth discussion on its applicability limits. In these 


works, the ID model was considered as a dimensional reduction of the 3D AE problem 
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near an isolated resonance. Recently, it was also demonstrated that holes and clumps 
may be generated with any (small) amount of background dissipation, provided that a 
phase-space plateau is formed by phase mixing and dissipative damping of an unstable 
kinetic resonance ( Lilley fc Nyqvist|2014 1. More precisely, in this case, holes and clumps 
are negative energy waves that grow because of background dissipation. The qualitative 
scenario of onset of stochastic transport in the presence of multiple resonances within 
the ID bump-on-tail paradigm has been recently reviewed by ([Breizman fc Sharapov 


2011 Breizman 20111 and the implications of quasi-linear diffusion in the presence of 


many modes have been discussed by (Lilley & Breizman 2012). In that work, it has 
been argued that extended flattening of particle distribution in the phase space can be 
more important than quasi-linear diffusion in the presence of many modes. The main 
limitation of these analyses is the adiabatic ansatz for nonlinear dynamics, i.e., the 
assumption that the characteristic nonlinear time is much greater than not 
only for quasi-linear diffusion but also for the evolution of chirping holes and clumps. 
Removing this limitation is one of the main motivation of the analysis proposed in the 
present work. Meanwhile, the conjecture that long-range frequency sweeping events can 
be mapped to corresponding AE behaviors in fusion plasmas is questionable, due to the 
lack of mode structures in the theoretical description ( Chen &: Zonca|[2015 ). Finally, we 
note that the Vlasov-Poisson system (on which the bump-on-tail paradigm is based) has 
been recently studied by numerical simulations based on Eulerian ( Shoucri|2010 | as well 
as Lagrangian codes (Pechhacker fc Tsiklauri|[2014|. 


3. Simulations of the two resonant beam model 


Let us now specialize the general problem described above to two supra-thermal beams. 
We first discuss the case when only the corresponding two most unstable waves are 
included in the dynamics; and we then address the problem of the presence of additional 
quasi-degenerate modes. It is worth noting that a similar problem was addressed by 


(Zaslavsky et al. 20081 in the framework of a Hamiltonian approach as in (Elskens & 


Escande 2003). However, that work includes an external magnetic field in the beam- 
plasma dynamics; resulting in a 3D geometry of the problem. Since the dispersion 
relation of the considered case contains additional degrees of freedom {e.g., transverse 
wave numbers), a direct comparison with the present analysis is limited, especially in the 
respect of considering two distinct resonant velocities as below. 


3.1. Most unstable modes 

The two most unstable modes are characterized by the resonance conditions ki = ojpfvi 
and ^2 = ujplv 2 - In what follows, we define 0 = V 2 /V 1 (thus implying £ 26 * = £ 1 ) and, 
since the system is specular with respect to the value (9=1, for the sake of simplicity, 
we report here only the results for 0 < 1 . In order to properly represent the phase-space 
dynamics, we move to the average beam speed reference frame. The dynamics of the two 
beam sheets is thus represented as 

Xu = Vit + Cli = t{vi + U 2)/2 -I- ^u{t) , X2z = V2t + C,2i = t{vi F U 2)/2 -£ ^2i{t) , ( 3 . 1 ) 

where it worth noting that the time evolution is formally the same for the different 
representations of beam nonlinear shifts but with modified initial conditions. Introducing 
the frequency mismatch as 


j3j = [fcj (ui -b U 2)/2 - ujp] /wp-q , 


(3.2) 
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Figure 1. Plots of the temporal evolution of \(j)j\ determined by the system ( |3.3[ ). We have set 
two different values of the parameter O as indicated in the graphs and the simulations are run 
for CTi = (72. The initial conditions and the parameter values are discussed in the text. 
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Figure 2. Reduced phase-space snapshots for fixed ti mes (indicate in the graphs). Each 

point corresponds to a charge sheet labeled by i in the system (3.31, the beam 1 (2) is represented 
in blue (red), the beams have equal densities. We have set & — 0.95 (left-hand panel) and 
0 = 0.97 (right-hand panel). (Color online) 


and using the scaled variables (2.10) (here ^ai = 27r^ai/L), the set of equations describing 
the nonlinear dynamics are 


= i (^-1 ^ ^-1 ^ 

= i (£-1 (j)^ ^ ^-1 ^ ^ 




L = ie-^^^^ 


2=1 

Ni 

(—T 

\Ni ^ 


No 


e 

N2 


E'■■'■'”) ■ 




(72 


N2 


,-it2?2. 


(3.3a) 

(3.36) 

(3.3c) 

(3.3d) 


In what follows, these equations are solved via a Runge-Kutta (fourth order) algorithm 
assuming -|- iV2 = 10^ total charge particle sheets and fj = 0.01. For the considered 
time scales and for an integration step h = 0.01, the total momentum and energy (derived 


from Eqs.(2.13l and (2.14) written in the average velocity reference frame) are conserved 
with relative fluctuations of about 1 x 10“^ and 4 x 10“®, respectively. In FIG{^ we plot 
the resulting wave amplitude, obtained as a function of time, setting (Ji = a 2 = 0.5 in 
order give equal weight to the two beams. As initial conditions, we consider the helds 
of order 0(10“^) and the sheets randomly distributed in 0 < ■Ci,2(0) < 27r, while the 
velocities at r = 0 ar e (0) = —^2(0) = — V 2 )/Luipfj = (1 — 0)j2£if]. Finally, the 

values of f)j from Eq.(3.2) are now f3j = [£j(l -|- Q)l2l\ — l]/^, with £2 = £\IQ- 

Erom the behavior of the mode amplitudes, the presence of a threshold value of 0 
clearly emerges, below which the two beams evolve independently generating two well 
separated single mode prohles. The transition to the interaction regime is anyway rather 
smooth. 

When 0 < 0.96, the generated field profile corresponds to two independent resonant 
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Figure 3. Plots of the temporal evolution of by setting 0 = 0.96 and cri = 0.8, a 2 = 0.2. 


modes. The amplification of both the waves proceeds as in the O’Neil analysis, where 
the initial exponential instability is followed by periodic oscillations. Furthermore, the 
instability is proportional to the corresponding a values consistent with linear theory. 
The dynamical independence of the beams can be argued also from the left-hand panel of 
FIGj^ where the phase space is depicted at a fixed time. Coherent with the single wave 
model, the dynamics of the system leads to the formation of rotating clumps which follow 
the corresponding Langmuir mode. In this sense, the number of such coherent structures 
is related to the wave number and we observe £i (beam 1 in blue) and £2 (beam 2 in red) 
clumps in the periodicity length L — 2n. The synchronized clumps evolve independently 
and no interference or transport between the two beams occurs. This behavior has been 
obtained in all the simulations for 0 < 0.96. 

When 0 > 0.96, the onset of interaction phenomena between the two beams is clear 
from FIG{^ Analyzing the phase-space dynamics in this case (in the right-hand panel of 
FIG[^ we set 0 — 0.97), we observe how during the formation of the rotating clumps, i.e., 
during the initial exponential evolution of the mode, the charge sheets of the first beam 
are involved in the dynamics associated to the second one. This yields the evolution of 
the Langmuir waves described in FIG for the same parameter values, where the leading 
mode is indeed (f> 2 - Such a mechanism is consistent with the fact that the mode (j )2 is 
excited by the slower beam and can tap energy from the mode (pi via wave-particle 
interactions with the faster beam particles. 

Interaction phenomena are evident also in the case of different values of the beam 
densities. In this respect, it is worth noting that a peculiar feature emerges near the 
threshold 0-value when cti is significantly greater than CT2, see FIG|| The larger density of 
the first beam forces the second mode to evolve differently compared to the equal density 
case. In particular, the weaker mode (p 2 is increased by the presence of the other one, 
again, because of energy transfer mediated by wave-particle interactions. This mechanism 
is analogous but opposite to the “surf” phenomenon occurring in plasma acceleration 
processes ( Litos et a/.|2014 Esarey et aZ.|1996 Keinigs fc Jones|1987 ). Specifically, when 
the second beam velocity V 2 is close to the enhanced Langmuir mode phase velocity, 
p 2 continues to tap energy from the faster beam. Thus, the weaker/slower mode can 
grow at the expense of the stronger/faster mode, via power transfer mediated by wave- 
particle interactions. In order to clarify the mechanism underlying this phenomenon, we 
can plot the power exchanged between each beam and individual modes, in the case 
under consideration (0 = 0.96 and cti = 0.8, CT 2 = 0.2). For the beam a and the mode j 
the dimensionless power can be written, using the scaled variables in the reference frame 
of the average beam speed, as: 


p(j) — 


Na 

E 

2 = 1 


O' 


)-l 


Jij£ai+iPjT 




7r(ui -I- V 2 ) 

uipTfL 


)+c.o.). 


(3.4) 
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Figure 4. Temporal evolution of (in the left-hand panels, the mode j = 1; in the right-hand 
panels, the mode j = 2) by setting 0 — 0.96 and cti = 0.8, (J 2 = 0.2. The power exchanged by 
the beam a = 1 {2) is plotted in black (red). In the upper (lower) panels, we have time averaged 
the power over a step At = 2 {At = 4) in order to smooth out small rapid oscillations. The 
transferred energy to the mode I 2 from the two beams is comparable, indicating the extraction 
of power from the mode j = 1 to the mode j = 2. (Colors online) 



Figure 5. Plots of the temporal evolution of \4>j\ by setting 0 — 0.998 and cti = 0.5, 02 = 0.5. 


This quantity is plotted in FIG|^ as a function of r separately for the two modes and 
with two different time averages as indicated in the figure, in order to better illustrate 
the global power transfer features. It is easy to realize from the right-hand panels that, 
when the “anti-surf” mechanism begins (r ~ 20, see FIGj^ the power exchanged by the 
two beams is comparable. Indeed, the particles of the fastest beam transfer energy to the 
mode £2 and drain energy from the mode £1. This behavior of power transfer confirms 
and clarifies the time evolution of the mode amplitude sketched in FIG|^ 

We conclude by briefly discussing the limit case for O of order fj. In particular, when 
& > 0.995 we observe the onset of a very particular behavior of the mode nonlinear 
evolution, as indicated in FIG[5]where we have set, as an example, O = 0.998. This is due 
to the fact that the beam instabilities overlap and the resulting dynamics is completely 
different from the previous cases. As soon as the clumps form, the beam charge-sheet 
interference is strong and the two beams can be effectively treated as a single phase-space 
structure but in the presence of peculiar time-dependent potentials generated by a fully 
interacting regime. 


3.2. Quasi-degenerate modes 

Let us now analyze how, in addition to the two most unstable waves discussed above, 
the two beam system can excite other nearly degenerate modes, i.e., fluctuations at 
frequency uij ~ tOp implying ~ 0. 
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Figure 6. First three panels from the left: Plots of the temporal evolution of \4>2\, 
and \4>s\ for different values of 0 as indicated in the graphs. Right panel: Plot of the mode (f>b 
as function of time. For this figure, we set cri = 02 and fj — 0.01. (Color online) 


For O over the threshold value for efficient nonlinear beam coupling, these modes 
correspond to sidebands, defined as = k 2 + pAk and kj = ki — pAk (where, Ak = 
k 2 — ki and p is an integer number). These quasi-degenerate modes are responsible for 
spectral broadening and, possibly, for a cascade phenomenon. To illustrate this process, 
we now discuss the case p = 1, considering thus the two sideband modes fcj" = 2 k 2 — ki 
and k^ = 2ki — k 2 - In the reference frame of the average beam speed, the dynamical 
equations governing the system can be easily constructed generalizing Eqs.(3.3) for the 
two additional quasi-degenerate modes. The corresponding mismatch parameters are 

= [i/e - 0/2 - 1 / 2 ] /^, ( 33 - = [0 - 1/20 - 1 / 2 ] /f). 

From the first three panels of Ficj^ only k'^ mode can be excited, while k^ does not 
grow since its resonance condition corresponds to velocities larger than vi and, further, 
it is associated with a negative value of the dielectric function (its frequency is actually 
slightly smaller than ujp). It is important to note that the unstable mode at k'^ saturates 
at comparable amplitude to the resonant mode ^2 although saturation is reached at 
a later time. Meanwhile, k'^ sideband excitation can occur only above the threshold 
condition in 0 for efficient nonlinear beam coupling. In fact, the mechanism for sideband 
excitation is the same as that for power transfer from faster to slower resonant modes 
discussed above. 

Efficient sideband generation can also occur in the case of two independent resonant 
beam modes, under special limiting conditions. This is the case when 0 —>■ 1/2“; i.e., 
k 2 —> 2k^ and the beat mode ki, = k 2 — ki —>■ k^. In fact, in this case, illustrated in 
the right-hand panel of EIGj^ the fcf, mode {4>b in the plot) is quasi-degenerate with 
the Langmuir frequency uip and has a resonance condition at velocities slightly below vi. 
Thus, it must be treated nonlinearly on the same footing of the two most unstable modes. 
As shown in the figure, the saturation level of the quasi-degenerate mode approaches the 
same value of 4>2 saturation, while the evolution of (pi is significantly altered because pb 
absorbs a relevant part of its energy via nonlinear wave-particle interaction. It can be 
also demonstrated that pb is not efficiently exited as 0 —?► 1. This fact is consistent with 
the original problem studied by ( O’Neil et aZ.||197T |, where mode-mode interactions are 
negligible. 


4. Lagrangian Coherent Structures and the transport processes 

In this Section, we analyze the transport processes in order to identify the regions 
where mixing phenomena are fastest and, in particular, their boundaries that are known 
as LCSs. The LCS technique was introduced by G. Haller in order to generalize intrinsi¬ 
cally asymptotic concepts in the dynamical system theory, such as stable and unstable 
manifolds, to the study of systems over a finite time span At. In the same way the 
separatrix of a pendulum splits the phase space into two macro regions where initial 
conditions have different asymptotic behavior, i.e., rotations and librations, the LGSs of 
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the beam-plasma system split the particles into two groups, i.e., recirculating and non¬ 
recirculating ones over the finite time span At. Without entering into the mathematical 


details of the definition, following Haller (Haller 20151 we state that LCSs need to be 
defined as material lines, i.e., lines advected by the flow, which are locally most repulsive 
with respect to the nearby ones over the time span At. Two initial conditions taken 
on different sides of an LCS will increase their distance over the time span At more 
than two initial conditions on two sides of any other curve. Therefore they will have, 
in general, a qualitatively different evolution the same way that two initial conditions 
starting on different sides of the stable manifold of a pendulum will increase their distance 
more than two initial conditions taken on the same side. This leads to the dehnition of 
a functional, called repulsion rate, which needs to be maximized on the LCS trough 


a variational approach (Haller 2011 Falessi et al. 20151. This technique has already 


been used in the computational fluid dynamics community ( 

Haller 

2001 

Shadden et al. 

2005 

Senatore & Ross||2011 Pascucci et aZ.||2010 

Garth et aL||2007 

) and, more recently, 

in plasma physics ( 

Borgogno et al. 

2011a|6). An introductory discussion on the proper 


operational definition of these structures can be found in (Falessi et al. 20151. Meanwhile, 


for a different approach to the characterization of coherent structures, see (Leoncini & 
Zaslavsk}^|2002 ). 


4.1. General approach 


In the recent review (Haller 2015), the author states that an LCS can be sought as 


an extremizing curve of the FTLE held, i.e., a ridge. Actually, a debate exists in the 


literature about the proper dehnition of a ridge. We refer the interested reader to (Haller 


2015) or to (Falessi et al. 2015) for more in-depth discussions of applications of this 


method in plasma physics. In our analysis, following ( Haller] 20151, we dehne a ridge of 
the FTLE held as a curve along which one sees a smaller change in the value of the held 
than in the direction transverse to it. The resultant contour plot of the FTLE is very 
peaked and almost constant on nearly ID regions. We will identify these structures with 


the ridges of the FTLE field and, therefore, with the LCSs of the system ( 

Lekien & Ross 

2010 Kasten et aZ.||2010 

Voth et oZ.||2002 

Chian et aL||2014 

Tang et a/.||2010 

I- 


In order to characterize the FTLE helds for a given conhguration, the trajectories of 
several test particles are traced when evolved in the time-dependent potentials r) 
(using the notation of the two-beam problem) generated self-consistently by the charge 
distribution. For a single reference time, the test particles are initialized to sample the 
whole phase space of interest (depending on the parameters). The system is composed of 
an arbitrarily large number of particles and, therefore, we can neglect the contribution of 
a single test-particle in the modihcation of the potentials, which are treated as assigned 
for a correspondent beam-plasma simulation. This approximation has already been used 
in ( Smith fc Pereira]|1978[ Tennyson et aL|[l994 |. 


To evaluate the FTLE profile for a hxed time r, the test particles are initialized in two 
phase-space grids of M/2 particles each and having an initial infinitesimal displacement 
in the velocity direction. In other words, a test particle located in has a neighbor 

initialized in [^, -f 6t]. With the assigned time-dependent potentials, at a time r -I- At 
the considered two test particles will be at a distance Sat in the phase-space and the 
FTLE value a in the point [^, f'j is evaluated as ( HallerpOll | 

(4.1) 


T, At) = In {5ATl5r)IAT . 


When considering a positive time evolution At > 0, the curves where the FTLE field is 
peaked define a repulsive transport barrier, while setting At < 0 and evolving the system 
backward in time they represent an attractive barrier. Approximated ID structures which 
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Figure 7. Overlay of the phase-space snapshots (black or red points) and the LCSs (thick lines) 
determine d as the peaks of the FTLE fields extracted from a corresponding contour plot built 
using Eq.(4.1|. Left panel: LCSs for the single wave model for r = 12.5 and At = 4. Middle 
panel: LCSsfor r = 11 and At = 10 (in this plot, the beam particles are highlighted in red). 
Right panel: LCSs in the late evolution for r = 156 and At = 4. (Color online) 


correspond to the LCSs can be built by plotting the maximum values of At) 

as extracted from a contour plot. In the analysis of this Section, we have set two grids 
of 250 X 250 test particles (M = 125000), thus obtaining 62500 values of the FTLE for 
each phase-space snapshot. The self-consistent fields ^(t) have been sampled from the 
complete Hamiltonian simulation with the proper selected parameters. 

4.2. The single wave model 

As first step, we apply the FTLE method to characterize the dynamical system 


describing the single wave model (O’Neil et aZ.||1971|, which can be easily deduced from 


Eqs.(2.12) by considering one beam and one mode j = 1. In (Tennyson et al. 1994| 


and (Smith & Pereira 19781, the test-particle approximation is used for studying the 


behavior of the single wave model using the dynamical systems techniques. In particular, 
in the latter, the evolution of the electromagnetic field is modeled with an analytic 
function while, in the former, it is extracted from the numerical simulations as in the 
present paper. In (Tennyson et al. 1994), the periodicity of the asymptotic state of the 


system is exploited in order to generate a Poincare section (stroboscopic map) of the 
motion, i.e., the sampling of the evolution of a particle after each periodic oscillation 
of the electromagnetic field. The set of points obtained by this procedure, for different 
initial conditions, can be used to distinguish different kinds of motion: periodic, quasi- 
periodic or chaotic. In both works mentioned above, a region of non-chaotic motion is 
present. This region moves, with a periodic behavior, coherently, i.e., only with slight 
deformations in its shape. Some of the beam particles will start inside this structure and, 
thus, they will be trapped during their entire evolution. This is the mechanism responsible 
for the formation of the clump in the Hamiltonian system and for the bunching of the 
particles. This procedure splits the ensemble of trajectories and, thus, the ensemble of 
the available particle motions, relying on their asymptotic behavior. For this reason, this 
scheme identifies barriers for the transport processes in the system which holds at any 
time, i.e., curves which cannot be crossed by the particles during their motion. In the 
most general aperiodic case however, such structures do not exist because every trajectory 
could in principle be chaotic and, furthermore, it is not clear how the stroboscopic map 
technique can be used for studying the single wave model during the linear growth. 

In this work, we are interested in studying the system during the onset of the instability. 
For this reason, the periodicity of the electromagnetic potential is generally lost. Using the 
FTLE method in order to identify the shapes of the LCS can nonetheless split the phase 
space into different regions based on the behavior of the trajectories over a finite period 
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of time. These structures, which are characterized by a peaked FTLE field, approximate 
the unstable (stable) manifolds of the dynamical system and characterize the transport 
through the lobes dynamics ( Malhotra fc Wiggins||1998 1 . 

The two parameters that characterize the FTLE values are r and At. The time 
evolution of the LCS is identified by the first, while the second determines how refined 
the LCS will be. By choosing an increasingly longer integration time At, the obtained 
LCSs acquire more intricate structures because of an increasingly better approximation to 


the homoclinic tangle. This kind of dynamics have been studied by (Malhotra & Wiggins 


19981 and it leads to the splitting of the phase space into macro-regions which exchange a 


relatively small amount of particles over a characteristic time scale. The characterization 
of the detailed shape of the stable and unstable manifolds, which can be achieved trough 
an FTLE plot with long At, gives insights into the time asymptotic transport processes. 
The shape of the LCS obtained with smaller values of At, meanwhile, highlights the 
boundaries of macro-regions with small transport among them. 

Here, we are interested in the shape and the location of these regions and, thus, the 
integration time is fixed as At = 4 after a number of tests, obtaining a simple, closed 
region as shown in the left-hand panel of FIGj^ where we have over-plotted the beam- 
particle positions to the LCS. The effect discussed above of increasing At can be deduced 
by comparing with the middle panel. Increasing the value of t, we get the position of 
the LCS at different instants and, therefore, the motion of the “finite time” transport 
barriers. From these plots, it can be argued that, during the onset of the instability, most 
of the beam particles are trapped inside the central region delimited by the LCS and, 
therefore, they recirculate and start forming the clump. 

In order to make a comparison with ( [Tennyson et al. 1994), we can calculate the LCS for 
larger values of r, when the electromagnetic potential oscillates periodically. We obtain 
the plot on the right-hand panel of FIGj^ where the transport barriers are splitting 
the chaotic region into two different areas having qualitatively different motions on a 
short time scale. The points outside the central “eye” have a motion similar to rotations 
in a pendulum, while the points inside the “eye” resemble librations. In contrast to an 
unperturbed pendulum, the points may start undergoing rotational motion and later 
undergo librational motion and vice-versa. However, this transition has a typical “long” 
time scale and, if we look at the transport processes which happen over a shorter time, the 
two different kinds of trajectories do not mix. For this reason the LCS may be regarded 
as finite time transport barriers. 


4.3. The two-wave model 

Let us now analyze the transport processes in the two-beam system using the same 
methodology introduced above. In FIGj^ we plot the LCS and the particle phase space 
at different values of O during the onset of the instability for r = 11.6 and At = 5. 

Considering the different colored particles initially belonging to different beams, it can 
be argued that they form a single phase-space structure for sufficiently small values of 
1 — 0. In the first plot, two eye-shaped regions delimited by the LCS clearly emerge. 
Correspondingly, the particles of the different beams do not mix. The shape of the LCS 
becomes increasingly complicated as 0 —>■ 1. In the most general case, the LCSs do not 
divide the phase space into separate regions and, correspondingly, the particles of different 
beams can mix into a single structure. Therefore, the mixing process is greatly enhanced 
when LCSs belonging to different beams intersect; and they cannot be considered as 
finite time transport barriers between separate phase-space regions. We refer to Ocrit as 
the smallest value of 0 compatible with the merging of the beams and give an estimate 
of such a threshold by increasing 0 until two LCS belonging to different beams intersect. 
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Figure 8. Overlay of the phase-space snapshots (blue and red points for the beam 1 and 
2, respectively) and th e LC Ss (thick lines) defined trough the contour plot peaks of the FTLE 
fields evaluated by Eq.( |4.1| |. The beam densities are set as cri = CT 2 = 0.5 and we gave considered 
T = 11.6 and At = 5. The different values of the beam initial velocity ratio (0) are indicated 
in the graphs. (Color online) 


From the plots in FIGj^ we obtain Ocrit between 0.96 and 0.97. For different beam 
densities, this method yields similar values of Ocrit- In the work ( Borgogno et aZ.|[2008 ), 
a similar analysis was carried out in order to study the transition from local to global 
chaos in the magnetic field lines of a reconnecting current layer. 

We conclude this Section with a derivation of Ocrit using the Chirikov criterion. 
Following (O’Neil et al. 1971), the velocity spread of particles trapped by a single 
saturated beam mode is Avj = V- The Chirikov parameter S is, thus, defined 

by the ratio between the resonance width and the resonance separation of the two beams 
with initial velocity vi and V2 {i.e., {vi — V2)/vi). It can then be written as 


2Avj _ 

(vi - V2)/Vi l-O 


(4.2) 


By the Chirikov criterion, two resonances overlap if S' > 1; therefore, in the saturation 
regime [i.e., (I)j ~ 1), we obtain Ocrit ~ 0.97, consistent with the previous analysis based 
on the LCSs (we remind that fj is fixed as 0.01). We note that, even if the definition 
of LCSs we have used is not rigorous nor unique, their shape is a posteriori consistent 
with particle dynamics and with the ensuing transport processes. The main advantage 
of using the LCS technique in studying transport processes between multiple resonant 
beams is that they can be calculated when the system is not near the saturated regime; 
thus, also when the Chirikov criterion cannot be straightforwardly applied. 


5. Concluding remarks 

In this work, we have proposed to investigate the ID bump-on-tail problem in terms of n 
cold beams injected in a cold background plasma behaving as a linear dielectric medium; 
and m ^ n self-consistently coupled nonlinear oscillators, mimicking the features of 
a warm fast particle distribution interacting with a spectrum of resonantly excited 
Langmuir waves. The use of the Hamiltonian formulation addressed in this work provides 
a convenient representation of the beam-plasma interaction and of ensuing transport 
processes at the disparate time and velocity scales underlying the complex behaviors of 
this system. 

Specializing the model to two supra-thermal beams, we have determined the threshold 
ratio between the initial beam velocities, for which the two corresponding instabilities 
mutually interact instead of remaining separate growing modes. Such a value has been 
found to be ~ 0.96 for the present simulation parameters, and is consistent with the 
estimates provided in (Levin et al. 1972) for a supra-thermal hot beam as well as with 
the Chirikov criterion. Above the threshold, the phase-space analysis of the particle 
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motion clearly shows that transport appears as an exchange of particles between the two 
beams. 

Furthermore, we have discussed the role of other quasi-degenerate modes, with partic¬ 
ular attention to sideband modes. More specifically, we have shown that such additional 
modes may nonlinearly saturate at almost the same level reached by the most unstable 
waves. This fact suggests that the long time-scale evolution of the multi-beam system 
would be significantly affected by nonlinear sideband generation; and the same would 
hold for EP transport. 

In the last par of the paper, we have investigated the onset of transport in the velocity 
space in terms of finite time transport barriers defined by the LCSs. The critical threshold 
for transport onset has been analyzed in terms of the LCSs of the system defined trough 
the FTLE field. The resulting critical value for the beam initial velocity ratio coincides 
with the Chirikov criterion and the direct inspection of particle motion in phase space; 
providing a very meaningful test of validity for the LCS technique in characterizing the 
phase-space dynamics. 

Our analysis provides a qualitative and quantitative investigation framework for phase- 
space transport processes in the beam-plasma system. The advantage in using the 
present paradigm stands in the direct and very general characterization of the phase- 
space morphology. This leads to the definition of proper characteristic time and velocity 
scales related with spectral features and wave-particle resonances, together with the 
characterization of the diffusion coefficient in the quasi-linear limit expected for a large 
number of modes and a sufficiently broad spectrum. 
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